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We study the role of finiteness and fluctuations about average quantities for basic structural 
properties of growing networks. We first determine the exact degree distribution of finite networks 
by generating function approaches. The resulting distributions exhibit an unusual finite-size scaling 
behavior and they are also sensitive to the initial conditions. We argue that fluctuations in the 
number of nodes of degree k become Gaussian for fixed degree as the size of the network diverges. 
We also characterize the fluctuations between different realizations of the network in terms of higher 
moments of the degree distribution. 

PACS numbers: 02.50.Cw, 05.40.-a, 05.50.+q, 87.18.Sn 



I. INTRODUCTION 



Networks such as the Internet and the World-Wide 
Web do not grow in an orderly manner. For example, 
the Web is created by the uncoordinated effort of millions 
of users and thus lacks an engineered architecture. Al- 
though such networks are complex in structure , their 
large size is a simplifying feature, and for infinitely large 
networks the rate equation approach H provides analyt- 
ical predictions for basic network characteristics. Never- 
theless, social and technological networks are not large 
in a thermodynamic sense {e.g., the number of molecules 
in a glass of water vastly exceeds the number of routers 
in the Internet). Thus fluctuations in network properties 
can be expected to play a more prominent role than in 
thermodynamic systems Qj . Additionally, extreme prop- 
erties, such as the degree the node with the most links 
in a network HH, the website with the most hyperlinks, 
or the wealth of the richest person in a society, are im- 
portant characteristics of finite systems. The size depen- 
dence of these properties or their distribution is difficult 
to treat within a rate equation approach. 

In this paper, we examine the role of finiteness and the 
nature of fluctuations about mean values for large, but 
finite growing networks. We shall focus primarily on the 
degree distribution iVfc ( N ) , the number of nodes that are 
linked to k other nodes in a network of N links, as well 
as related local structural characteristics. We shall argue 
that self averaging holds for the degree distribution, so 
that the random variables Nk(N) become sharply peaked 
about their average values in the N — > oo limit. We 
shall also argue that the probability distribution for the 
number of nodes of fixed degree, P(Nk,N), is generally 
a Gaussian, with fluctuations that vanish as N — ► oo. 
On the other hand, higher moments of the degree dis- 
tribution do not self average. This loss of self-averaging 
ultimately stems from the power-law tail in the degree 
distribution itself. 

In the next section, we define the growing network 
model and briefly review the behavior of the average de- 
gree distribution in the thermodynamic N — > oo limit. 
We also discuss how the average degree distribution can 
naturally be expected to attain a finite-size scaling form 



for large but finite N. We then describe our general 
strategy for studying fluctuations in these growing net- 
works. In Sec. Ill, we outline our simulational approach 
and present data for the average degree distribution. In 
the following two sections, we examine the role of finite- 
ness on the degree distribution, both within a continuous 
formulation based on the rate equations (Sec. IV), and an 
exact discrete approach (Sec. V). The former approach is 
the one that is conventionally applied to study the kinet- 
ics of evolving sustems, such as growing networks. While 
this approach has the advantage of simplicity and it pro- 
vides an accurate description for the degree distribution 
in an appropriate degree range, it is quantitatively inac- 
curate in the large degree limit. This is the domain where 
discreteness effects play an important role and the exact 
discrete recursion relations for the evolution of the degree 
distribution are needed to fully account it properties. In 
Sec. VI, we discuss the implications of our results for 
higher moments of the degree distribution and their as- 
sociated fluctuations. Sec. VII provides conclusions and 
some perspectives. Calculational details are given in the 
appendices. 



II. STATEMENT OF THE PROBLEM 



The growing networks considered in this work are built 
by adding nodes to the network one at a time according 
to the rule that each new node attaches to a single pre- 
vious node with a rate proportional to A^, where k is 
the degree of the target node. We investigate the class of 
models in which Ak = k + A, where A > — 1, but is oth- 
erwise arbitrary. The general situation of — 1 < A < oo, 
corresponds to linear preferential attachment, but with 
an additive shift A in the rate. This model was origi- 
nally introduced by Simon to account for the word fre- 
quency distribution 0. The case A = corresponds to 
the Barabasi- Albert model [||, while the limit A — > oo 
corresponds to random attachment in which each node 
has an equal probability of attracting a connection from 
the new node. Thus by varying A, we can tune the rela- 
tive importance of popularity in the attachment rate. 

Previous work on the structure of such networks was 
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primarily concerned with the configuration-averaged de- 
gree distribution (N k (N)), where the angle brackets de- 
note an average over all realizations of the growth pro- 
cess for an ensemble of networks with the same initial 
condition. Additionally, most studies focused on the tail 
region where k is much smaller than any other scale in 
the system. For attachment rate A k = fc + A, this average 
degree distribution has a power-law tail 



(N k (N)) = Nn k , with n fe oc£:-( 3 + A ) 



(1) 



as iV — ► oo. In the specific case of A k — fc, the average 
degree distribution explicitly is |t] Q 



(N k (N)) = Nn k , with n k 



k(k + l)(k + 2) 



(2) 



For finite N, however, the degree distribution must 
eventually deviate from these predictions because the 
maximal degree cannot exceed N. To establish the range 
of applicability of Eqs. ([[]), we estimate the magnitude of 
the largest degree in the network, /c max by the extreme 
statistics criterion Y^, k > k (Nk{N)) « 1. This yields 
fcmax oc A^ 1 ^ 2+A ^. We therefore anticipate that the aver- 
age degree distribution will deviate from Eq. (|l]) when k 
becomes of the order of fc max . The existence of a maximal 
degree also suggests that the average degree distribution 
should attain a finite-size scaling form 



(N k (N)) ~ Nn k F(0, £ = k/k n 



(3) 



Some aspects of these finite-size corrections were recently 
studied in Refs. |L2] -|l5). One basic result of our work is 
that we can compute the scaling function explicitly. We 
find that this function is peaked for k of the order of fc max 
and that it depends substantially on the initial condition. 
In contrast, the small-degree tail of the distribution - the 
reason why such networks were dubbed scale-free - is in- 
dependent of N and the initial condition. 

To study finite networks where fluctuations can be sig- 
nificant, we need a stochastic approach rather than a de- 
terministic rate equation formulation. For finite JV, the 
state of a network is generally characterized by the set 
N = {Nx,N 2 ,...} that occurs with probability P(N). 
The network state N evolves by the following processes: 



(Nx,N 2 ) 
(Ni,N k ,N k+1 ) 



(Nx,N2 + l) 
{Nx + l,N k 



l,N k+1 + l). 



The first process corresponds to the new node attaching 
to an existing node of degree 1; in this case, the number 
of nodes of degree 1 does not change while the number 
of nodes of degree 2 increases by 1. The second line ac- 
counts for the new node attaching to a node of degree 
k > 1. 

From these processes, it is straightforward, in princi- 
ple, to write the master equation for the joint proba- 
bility distribution -P(N). It turns out that correlation 



functions of a given order are coupled only to correla- 
tion functions of the same and lower orders. Thus we 
do not need to invoke factorization (as in kinetic theory) 
and we could, in principle, solve for correlation functions 
recursively. However, this would provide much more in- 
formation than is of practical interest. Typically we are 
interested in the degree distribution, or perhaps two- 
body correlations functions of the form (iVjiVj). Even 
though straightforward in principle, it is difficult to com- 
pute even the two-point correlation functions (iVj Nj) for 
general i and j. In this work, we shall restrict ourselves 
to the specific (and simpler) examples of (Nf ), (iVi N2), 
and (iV|). We will use these results to help characterize 
fluctuations in finite networks. 



III. SIMULATION METHOD AND DATA 

To simulate a network with attachment rate A k — k+X 
efficiently, we exploit an equivalence to the growing net- 
work with re-direction (GNR) ||. In the GNR, a newly- 
introduced node n selects an earlier "target" node x uni- 
formly. With probability 1 — r, a link from n to x is cre- 
ated. However, with probability r, the link is re-directed 
to the ancestor node y of node x (Fig. [I]). As discussed 
in || , the GNR is equivalent to a growing network with 
attachment rate A k = k + A, with A = r -1 — 2. Thus, for 
example, the GNR with r ~ 1/2 corresponds to the grow- 
ing network with linear preferential attachment, A k = k. 
Simulation of the GNR is extremely simple because the 
selection of the initial target node is purely random and 
the ensuing re-direction step is local. 




probability 1-r 



probability r 

FIG. 1. The re-direction process. The new node n selects 
a random target node x. With probability 1-ra link is es- 
tablished to this target node (dashed), while with probability 
r the link is established to y, the ancestor of x (solid). 

There is, however, an important subtlety about this 
equivalence that was not discussed previously in Ref. [^) . 
Namely, the redirection process does not apply when a 
node has no ancestor. By construction, every node that 
is added to the network does have a single ancestor, but 
some primordial nodes may have none. For example, for 
the very natural "dimer" initial condition a — o, the seed 
node on the left has no ancestor and the GNR construc- 
tion for this node is ambiguous. One way to resolve this 
dilemma is to adopt the "triangle" initial condition in 
which there are 3 nodes in a triangle with cyclic con- 
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nections between nodes. This leads to the correct at- 
tachment rate for each node for any value of A. We 
therefore typically use this initial state to generate de- 
gree distribution data. On the other hand, theoretical 
analysis is simpler for the dimer initial condition. This 
state can also be simulated in a simple manner (for the 
case A = 0) by a slightly modified GNR construction in 
which direct attachment to the seed node is not allowed. 
It is straightforward to check that this additional rule 
leads to the correct attachment rates for all the nodes in 
the network. 
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FIG. 2. Normalized degree distributions for the triangle 
initial condition for networks of 10 2 , 10 3 , . . . links (upper left 
to lower right), with 10 J realizations for each N, for (a) 
A k = k (up to 10 6 links) and (b) A k = k + A, with A = -0.9 
(up to 10 links). In (a), the dashed line is the asymptotic 
result rik = 4/[k(k + l)(fc + 2)]; the last three data sets were 
averaged over 3, 9, and 27 points, respectively. In (b), the 
last two data sets were averaged over 10 and 100 points, re- 
spectively. 

Figure ^| shows the average degree distribution for at- 
tachment rates Ak — k and Ak = k + \ with A = —0.9 for 
the triangle initial condition. This latter value of A gives 
results that are representative for values of A close to — 1 . 



The data exhibits a shoulder at k 
more pronounced when A < (Fig. ||(b)). This shoulder 
is also at odds with the natural expectation that the aver- 
age degree distribution should exhibit a monotonic cutoff 
when k becomes of the order of k max . This shoulder turns 
into a clearly-resolved peak that exhibits relatively good 
data collapse when the degree distribution is re-expressed 
in the scaling form of Eq. (||) (Fig. ||). Conversely, the 
magnitude of the peak diminishes rapidly when A is pos- 
itive and becomes imperceptible for A > 0.5. 
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FIG. 3. The corresponding scaling function F(£) in Eq. (^) 
for the data in Fig. 0(a). 

In the following two sections, we will attempt to under- 
stand this anomalous feature of the degree distribution 
by studying the rate equations for the node degrees of 
finite networks. 



IV. CONTINUUM FORMULATION 

We focus on the case of the linear attachment rate 
Ak = k and briefly quote corresponding results for other 
attachment rates. In the continuum approach, N is 
treated as continuously varying. Then the change in the 
average degree distribution satisfies the rate equation 



d(N k (N)) / (k - l}N k -i(N) 



dN 



2N 



(4) 



We assume the dimer initial state - two nodes connected 
by a single link so that Nk(N = 1) = 2dki- 

Equations (^J) are recursive and can be solved sequen- 
tially, starting with (Ni). Explicit results for (Nk), k < 4, 
are given in Appendix A. These expressions show that the 
dominant contribution in the N — > oo limit is linear in ./V 
and this corresponds to the solution in Eq. (g). Indeed, 
if we substitute (Nk{N)) — n k N into Eqs. (Q), we obtain 
the recursion nf. = rik-i{k — l)/(fc + 2), whose solution 
is Eq. (^). From the first few {Nk}, it is easy to see that 
the first correction to this leading behavior is of the order 
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of N- 1 ' 2 . Substituting (N k (N)) = n k N + A k N- 1 ? 2 into 
Eqs. and keeping the first two terms in each (Nk), we 
find A k = 4/3. Continuing this procedure systematically, 
we arrive at the expansion: 



(MN)} 



1 



nkN + tm/ 2 

4 (k-l)(k - 2) 
5 
5 



3 k 
2 ~ 



f 



N 3 / 2 
(fe — l)(fe — 2)(fe 



N 



3) 



18 



N 2 



1 (fc-l)(fc-2)(fc-3)(fc-4) 

+ li WT 2 



(5) 



In general, the right-hand side contains k + 1 terms which 
can be written more succinctly as 



(N k (N))=n k N- 



1 



fc-i 



r(fc) (-1)^ 



Nil* 



(6) 



The coefficients ia,- = (2j+4)/[j!(j + 3)] may be obtained 
by imposing the initial condition N k {\) = 28 k ,i as each 
(N k ) is computed; a simpler way of obtaining these co- 
efficients will be explained below. Note that expansion 
(||) is asymptotic because successive terms decrease only 
when k -C y/~N. 

A more convenient way to solve Eqs. (jj) is in terms of 
the generating function 

oo 

J^(N 1 z) = ^(N k (N))^ (7) 

Multiplying Eq. (Q) by z k and summing over k, the gen- 
crating function satisfies the following partial differential 
equation 



2iV 9iV +z(1_ '" ) ^ )^(iV,2)=2iVz. (si 



The initial condition is 7V(1, z) = 2z, corresponding to a 
starting point of two nodes and a single connecting link. 

We reduce Eq. (||) to a wave equation with con- 
stant coefficients by changing from the variables (N, z) 
to (In y/N, ln[z/(l — z)]). Then by introducing the ro- 
tated coordinates x,y such that x + y = ]ny/~N and 
x — y — ln[z/(l — z)], we recast the wave equation into 



dN(x, y) 2 e 3x+2y 



dx 

whose general solution is 



(9) 



Af(x, y) = e 2x+2y - 2e x+3y + 2e 4y In (e* + e y ) + G(y). 

Finally, G(y) is found by imposing the initial condition 
Af(l,z) = 2z. When N = 1, we have x — —y, so that 
the initial condition becomes Af(—y,y) = 2/(1 + e 2y ). 
Therefore 



G(y) = 



and finally 



1 + e 2y 



1 + 2e 2y - 2e 4y In (e~ v + e y ) 



1 - p 2y 

N{x,y)=e 2y (e 2 *-2e x + y + 2)+ l Y ^ ry 

/ P x+V i ply \ 

+ 2e^ln(^±^). (10) 

Using e 2y — (z^ 1 — 1)>/N and e x+v = y/~N, we re-express 
the generating function in term of the original variables 



Af(N, z) = (3 - 2z- x )N + 2(z~ 1 
1 - {z- 1 - i)Vn 



1)VN 



1 + (z- 1 - 1)VN 
2(z~ 1 -l) 2 N ln( l-z + -^=] 



(11) 



We are primarily interested in the degree distribution 
for nodes whose degree is of the order of fc max » y/~N. 
This part of the distribution can be extracted from the 
limiting behavior of the generating function Af(N, z) as 
z — > 1 from below. Since the interesting range is k m \/~N , 
it is convenient to write 



1 



N 



(12) 



and keep s finite while taking N — > oo limit. We simplify 
still further by eliminating the contribution to the gen- 
erating function from the power-law tail of n k in Eq. (|^). 
For this purpose we consider the modified generating 
function 



d 



Af = Y l (k + 2)(k + l)k(N k ) z 



fe+3 



(13) 



fc=l 



which is constructed so that the derivatives multiply the 
degree distribution by just the right factors to eliminate 
the power law tail. The leading behavior of this modi- 
fied generating function will therefore provide the scaling 
function F(£) of Eq. (§). 

We now substitute Eq. (|lj) and the anticipated scaling 
form of Eq. (|) into the right-hand side of Eq. (|l3]) and 
replace the sum by an integral. This gives the Laplace 
transform of the scaling function times a prefactor, 



47V 3/i 



(14) 



with £ = kjN 1 ! 2 . Using Eq. (fill), we compute the deriva- 
tive on the left-hand side of Eq. (|13]). In the N — > oo 
limit, this derivative becomes AN 3 ' T J(s) with 



J(e) = 



1 



1 



l + a (i + s) 2 (i + s) 3 (i + s y 



(15) 
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This is just the Laplace transform of the scaling function. 
Inverting the Laplace transform then yields 



F(t) = (l + 0(l + f ) < ; 



(16) 



Notice that the coefficients Vj in Eq. (|^) can be obtained 
by expanding F in a Taylor series. This is a much simpler 
approach than solving each (Nk(N)) directly. 

An important feature of the degree distribution is that 
it depends significantly on the initial condition. For ex- 
ample, for the triangle initial condition, solving Eq. (^) 
subject to N A (N = 3) = 34,2, or Af A (3,z) = 3z 2 , yields 

M A (N, z) = (3 - 2z- 1 )N + 2{z- x - 1)V3N 



+ 3(l + (z- 1 -l) v /iV/3)" 



2(z~ 



1) 2 N In 1-z + z 




Repeating the steps used to deduce the scaling function 
(O) from Eq. (pT|), we now find 



F*(£)= 1 



, v 2 , v A 



v 



tVz. (18) 



Therefore small differences in the initial condition trans- 
late to discrepancies of the order of y/N in the degree 
distribution of a finite network of N links. Thus the 
properties of the nodes with the largest degrees are quite 
sensitive to the first few growth steps of the network (see 
also Ref. §). 

While this initial condition dependence is real, there is 
also a spurious aspect to this effect. This may be illus- 
trated by considering the linear trimer initial condition 
o — o — o. This is the unique outcome of the dimer initial 
condition after one node has been added. These two ini- 
tial conditions should therefore lead to the same degree 
distribution. However, for the linear trimer initial state 
(N k (N = 2) = 26^1+Sf. 2) the continuum approach gives 
the scaling function, 



F(0 



1 



n 



4 



V 



which is distinct from Eq. (jig)! This anomaly highlights 
one basic limitation of the continuum formulation. 

Finally, we mention that parallel results can be ob- 
tained for the general case of the shifted linear attach- 
ment rate, A k = k + X. The rate equation for the average 
degree distribution is 



d(N k (N)) _ I A k _ 1 N k _i(N) — A k N k (N) \ 
dN \ A J 



+ 5 



k,l 



where A = Y, A k N k = £)(fc + X)N k . To compute A we 
use the sum rules kN k = 2N (every link contributes 
twice to the total degree), as well a,s N k = N + 1 (for 



any tree initial condition) or ^2 N k — N (for an initial 
condition that has the topology of a single cycle) . To sim- 
plify final formulae, we use the latter topology (specifi- 
cally, the triangle initial condition) so that A = (2 + X)N. 

Solving the above rate equations successively, we find 
that the first two terms in the asymptotic series for 
(N A (N)) are 



(7V fc A (A0) ~ n k N + n' k N -(l+»/V+» 



(19) 



with 



n k 



r(3 + 2A) r(fc + A) 



+ A) r(l + A) T(k + 3 + 2X) ' 
2 + A 3 ( 3+2A )/( 2+A ) T(k + X) 



k 3 + 2A L(l + A) T(fc) ' 

The corresponding leading behaviors are n k oc fc~( 3 + A ) 
and n' k oc k x . Thus the two contributions to the de- 
gree distribution in Eq. ( |l9| ) are comparable when k ~ 
^yi/(2+A)^ This value coincides with maximal degree 
&max that is obtained by the extreme value condition 
Sfe>fe N/k 3+x w 1. Once again the degree distribu- 
tion is described by a scaling function in the dimension- 
less variable £ = fc/7V 1 /( 2+A ). 



V. DISCRETE APPROACH 

We now turn now to the discrete approach for the net- 
work evolution. That is, one link is introduced at each 
discrete time step; this corresponds exactly to what oc- 
curs in the simulation. We again focus on the caes of 
the linear attachment rate A k — k. We first treat in de- 
tail the case of nodes of degree one and then extend our 
approach to nodes of higher degrees. Finally, we give a 
scaling description for the degree distribution itself. 



A. Nodes of Degree One 

The number of nodes of degree one, N±(N), is a ran- 
dom variable that changes according to 

Pr0b ' 2^ 

2N 

after each node addition event. That is, with probabil- 
ity Ni/2N, a newly-introduced node attaches to a node 
of degree one; in this case, the number of nodes of de- 
gree one does not change. Conversely, with probability 
(1 — N\/2N), the new node attaches to a node of degree 
greater than one and N\ thus increases by one. Therefore 




prob. 1 



(20) 



N?(N)\ 



2N 



7 



JVi(JV) + l 



N?{N) 
2N 



Y,(A0\ 



2N 



5 



from which 

(N 1 (N+l)) = l+(l-^(N 1 (N)). (21) 

We take the initial condition (Ni(l)) = Ni(l) = 2. 

We solve this recursion in terms of the generating func- 
tion X x {w) = T,n>i( N i( N )) wJV_1 - We therefore multi- 
ply Eq. ( fj"l| ) by Nw 1 ^^ 1 and sum over N > 1 to convert 
this recursion into the differential equation 



dX x _ 1 1 dXi 

dw (1 — w) 2 2 1 dw 



(22) 



Solving Eq. (|22| ) subject to the initial condition Xi(0) = 2 
gives 



Xx{w)= 2 - ' 



1 



(23) 



3 {1-wf 3 (1-w) 1 / 2 ' 
Finally, we expand X\ (w) in a Taylor series in w to obtain 



2 4 T (N - -) 



(24) 



The leading term is identical to that in the continuum 
approach (cf. Appendix A), but the coefficient of the cor- 
rection term is 4/(3-^/7?) ~ 0.7523, compared to 4/3 in the 
continuum approach. 

The discrete approach is also suited to analyzing higher 
moments of the random variable Ni (N) . The second mo- 
ment (Nf(N)) plays an especially important role as we 
can then obtain the variance a\ = (N 2 (N)) - (Aq(AT)) 2 
and thereby quantify fluctuations. From Eq. (^oj) this 
second moment (N 2 (N)} obeys the following recursion 
formula 



(N 2 (N + 1)) =1+^1-1^ (N? 



The solution to this recursion is outlined in Appendix B 
and the final result is 



4 r(iV-l) 35 
lU/n T{N) Y ' 1 ' 



(26) 



In the large N limit, we use Stirling's formula to give, for 
the variance, 



N 20 1 



16 1 



9 90? N 1 / 2 9y^K N 



(27) 



To obtain the entire probability distribution P(Ni, N) 
one must solve 



P(N U N+1) = ^P(N 1 ,N) 
Ni - 1 



1 



2A^ 



P(Ni-l,N). (28) 



By the Markov nature of the process, P(N\,N) should 
approach a Gaussian distribution in the large N limit. 
Numerically, we indeed find a Gaussian distribution with 
a peak at 27V/3 and dispersion i y/N in agreement with 
our theoretical results for (Ni(N)} and (N 2 (N)). 



B. Degree Greater Than One 



For k > 2, the random variable Nk = Nk(N) changes 
according to 



kN k 



N k (N + l) 



N k - 1 prob 
N k + l prob 
N k prob. 1- 



2N 

(fc - 1)^_! 

27V 

(k - l)N k -i + kN k 
2N 



(29) 



at each node addition event. Again, because of the 
Markov nature of this process, we anticipate that 
P(N k ,N) approaches a Gaussian distribution for every 
fixed degree fc; therefore, we only need calculate (N k (N)} 
and (N 2 (N)) to infer the asymptotic distribution. To de- 
termine the first moment, we repeat the steps described 
in detail for fc = 1 and obtain the recursion formula 



(N k (N + l)) 



(N k (N)) 
(fc - l)N k - 



2N 



kN k (N) \^ 



(30) 



The solution to this recursion is given in Appendix C 
and explicit formulae for (N k (N)) for fc < 5 are also 
quoted. Qualitatively, these results closely correspond 
to the asymptotic series for (N k (N)} in the continuum 
formulation (Eq. ^)) but with somewhat different coef- 
ficients in the correction terms. 

The determination of the second moment (-/V|) is more 
complicated because it is coupled to (Nk-xNk), which 
in turn is coupled to (Nk-2Nk), etc. However, we can 
still determine (N 2 ) for small fc (Appendix D). From the 
structure of the rate equations, our general conclusion is 
that u\ = (N 2 (N)) - (N k (N)) 2 = fi k N. Therefore the 
distribution of N k (N) approaches a Gaussian for each fc 
as N — + oo. 
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C. Generating Function Approach 



D. Scaling Function 



In close analogy with Sec. Ill, we now obtain the gener- 
ating function for (Nk(N)), from which the exact scaling 
function in Eq. (^) can be deduced. Since Eq. (^0|) in- 
volves two discrete variables, k and N, it proves useful 
to introduce the two-variable generating function 



w N - 1 z k 



(31) 



N=l fc=l 



The governing equation for Af(w, z) that is obtained from 
Eq. (H), is 



2(1 _ m) £ + 2(1 _ 2) |L_ 2 )^ = _^_. (32) 



This is similar to Eq. (||) and can be solved accordingly. 
We introduce the rotated variables x, y such that 



1 z 
x + y = -- ln(l - w), x-y = ln——, 



to recast Eq. (|32| ) into 
d 



o 5x+4y 



e» 



(33) 



(34) 



The general solution is 



Af(x, y) = e ix+iy - 2e 3x+5y 

+ 2e 2x+6y In (e x + e v ) + e 2x G(y), 

and the function G(y) is found from the initial condi- 
tion ]\[{w = 0, z) = 2z. When w — 0, we have x = 
-y = | ln[z/(l-z0], and hence Af(-y, y) = 2/(1 + e 2y ). 
Therefore 



G(y) = 



2e 2y 
1 + e 2y 



- 2y + 2e iy - 2e 6y In ( e -*+e»} 



and finally 



Af(x, y) = e 4x+4y - 2e 3x+5y - e 2x+2y + 2e 2x+ 

p 2x+2y / x+y , 2y 

+ 2- ^ + 2e 2x+6y ln' 



4 y 



l + e 2y ' \ 1 + e 2y 

In term of the original w, z variables, 
(3-2Z- 1 ) 1 



(35) 



Af(w, z) 



(1 — w) 2 1 — w 
2(z- 1 -l) 2{l-w)- 1 ' 2 



(l_ w )3/a ( z -i _ !) + (! _ w y/2 
2{z- 1 -l) 2 



(1-w) 2 



In 



1 — z + z(l — It)) 



1/2 



(36) 



By expanding N{w, z), we can in principle determine all 
the (N k (N)). 



To extract the scaling function F(£) from the gener- 
ating function N(w, z) we use the same approach as in 
Sec. IV. The details are given in Appendix E and the 
final result is 



= erfc 



2£ + e 3 



■in 



(37) 



where erfc(a;) is the complementary error function. A 
similar result for a related network model was found pre- 
viously by Dorogovtsev et al. [[uj. Notice that the exact 
form for F(£) vanishes much more quickly than predicted 
by the continuum approach. When k 3> y/N, the contin- 
uum approach gives 



(N k (N)) cont . -» — e- k /^ 



N 



(38) 



while the exact average degree distribution has a Gaus- 
sian large-degree tail 



'ttN 



-k 2 /4N 



(39) 



The scaling function in Eq. ( p7| ) quantitatively ac- 
counts for the shoulder in the degree distribution. In 
contrast, while the scaling function from the continuum 
approach does exhibit a peak, it is both quantitatively 
and qualitatively inaccurate (Fig. ^). 




k/N 



FIG. 4. Comparison between the scaling function F(£], 
with £ = k/N 1 / 2 , in the continuum approximation [Eq . (p_S|) , 
dashed curve] and in the discrete approach [Eq. (|37|), solid 
curve]. The circles give the simulation data of 10 realiza- 
tions of a network with N = 10 4 links for the dimer initial 
condition; these data coincide with the theoretical prediction. 
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VI. HIGHER MOMENTS AND THEIR 
FLUCTUATION 

We now turn to higher moments of the degree dis- 
tribution, as well as the fluctuation in these quantities 
between different realizations of the network. While the 
zcroth and first moments of the degree distribution are 
simply related to the total number of links for any net- 
work topology, the higher moments are not so simply 
characterized, but instead reflect the power-law tail of 
the degree distribution. 

We first compare the moments of the average degree 
distribution to appreciate the difference between the con- 
tinuum and exact descriptions. For the second moment, 
we use the identity 



£fc(fc+l)(iV fe 



k=l 



d_ 

dz 



Af(N, z) 



(40) 



Using N{N, z) from Eq. (|Tl]), together with the value of 
the first moment, we obtain, in the continuum approxi- 
mation, 



(k 2 



2NlnN + 2N. 



(41) 



fc=i 



On the other hand, using the exact discrete expression 
(pw we find 



d_ 

dz 



2 = 1 



- 21n(l - w) 
(1 - w) 2 



which we then expand in a series in w to yield, for the 
second moment, 



(k 2 



= 2NH 



N ■ 



(42) 



fe=i 



Here Hn = ^2 1< j <N j 1 is the harmonic number (jl6| 
In the large N limit, therefore, 



M, 



/(fc™)\ 
\(k°)J 



l/n 



(45) 



has the following TV dependence: 
M n oc 



const. n < 2 

In N n = 2 

N (n-2)/2 n>2 



(46) 



In a related vein, we also study the fluctuations in these 
moments between different realizations of the network 
growth. That is, we record the value of (k 2 ) for each 
realization of the network to obtain the underlying dis- 
tribution P((k 2 )). A typical result is shown in Fig. [jj. 
Notice that the distribution of (k 2 ) is relatively broad 
with an exponential tail. The distributions of higher mo- 
ments are even broader, with each being dominated by 
the realizations with the largest value of the correspond- 
ing moment. 



10" 



10" 



10" 



10 



10 



10000 



20000 

<k 2 > 



30000 



FIG. 5. Distribution of (k 2 ) for 10 5 realizations of a grow- 
ing network with N — 10 3 for attachment rate Ak = k 
with the triangle initial condition. The raw data has been 
smoothed over a 100-point range. 



(k 2 



2N\nN + 2jN + 1 



67V 



where 7 = 0.5772166 is Euler's constant. 

For higher moments, even the leading term given by 
the continuum approach is erroneous. For example, 

(fc 3 ) cont . = 24iV 3 / 2 - 67V In N - 22N , (43) 

while the exact value is 



(fc 3 



32 T(N- 



T(N) 



6NH N -WN. 



(44) 



More generally, the dependence of the moments on N 
stems from the power-law tail of the degree distribution 
(Nk) oc N/k 3 . From this asymptotic distribution, a suit- 
ably normalized set of measures for the mean degree 



VII. CONCLUDING REMARKS 



We studied the role of finiteness on the degree distribu- 
tions of growing networks with a node attachment rate of 
the form Ak — k + X. For finite networks, fluctuations are 
no longer negligible and a stochastic approach is needed 
to analyze these properties. We found the average degree 
distribution within an approximate continuum formula- 
tion and by an exact discrete approach. The continuum 
approach has the advantage of being much simpler than 
the discrete formulation, but does not provide a quan- 
titatively accurate description of the large-A; tail of the 
degree distribution. 

We also argued that the degree distribution Nk(N), 
when considered as the random variable in k, exhibits self 
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averaging, i.e., the relative fluctuations in Nk(N) dimin- 
ish as N — ► co. Moreover, the variance <j\ = (N^)-(Nk) 2 
varies linearly with N, and the probability distribution 
P(Nk,N) approaches a Gaussian. To support these as- 
sertions, we computed a 2 for k — 1,2. These partial 
results support our general hypothesis that fluctuations 
in Nk(N) are Gaussian. Perhaps the Van Kampen f2- 
expansion |fl7|| would prove to be a more appropriate 
analysis tool to undertake a systematic study of fluctua- 
tions in growing networks. 

Of course, the random variables JVfc(JV) should be 
Gaussian only for sufficiently small k, viz., as long as 
(N k ) > 1, or equivalently, k < iV 1 /(3+A)_ Qn the other 
hand, fluctuations become large and non-Gaussian when 
k oc N X ^ 3+X K Determining the fluctuations in this de- 
gree range seems to be difficult, as one must study the 
master equation for the joint probability distribution. 

In this work, we limited ourselves to the degree dis- 
tribution; this is perhaps the most important and also 
the most easily analyzable local structural characteristic 
of a network. However, recent investigations of growing 
networks has increasingly focused on global characteris- 
tics, such as the size distribution of connected compo- 
nents, see e.g., Refs. |p^-|2l|]. The methods described 
in this paper should be applicable to probing fluctua- 
tions of the component size distribution and other global 
network characteristics. This direction seems especially 
exciting since the simplest growing network models that 
allow for a multiplicity of clusters exhibit a very unusual 
infinite-order percolation transition |l^-|2l|. Thus one 
might anticipate interesting giant fluctuations near the 
percolation transition of these models. 

We are grateful to NSF grant DMR9978902 for partial 
financial support of this research. 



APPENDIX A: THE AVERAGE DEGREE 
DISTRIBUTION IN THE CONTINUUM 
FORMULATION 



Within the continuum framework, the average degree 
distribution is described by Eqs. (|J). Successively solv- 
ing these equations by elementary methods, we obtain 
(N k (N)). For k = 1,2,3, and 4 we obtain: 



1 4 1 

(N 2 (N)) = -N + 

v v ;/ 6 3 N 1 ' 2 



3 1 

2 N' 



(N 3 (N)) = ^N+i 
(N4(N)) = ^N 



3 ATi/2 N 5 TV 3 / 2 ' 



1 



3 iVV2 



9 1 

2 N 



24 1 

~5 N^ 2 



5 1 
3 TV 2 "' 



APPENDIX B: GENERATING FUNCTION FOR 

(N?(N)) 

To determine (N 2 (N)), we introduce the generating 
function y t (w) = J2n>i( N i ( N )) This converts 

the recursion relation Eq. ( p5|) into the differential equa- 
tion for the generating function 



dyi 



i 



(| -''^=CT^ + 2* + 2 ^' ™ 



with X\{w) given by Eq. (J23|) . Solving ( |Bl| ) subject to 
the initial condition 3^i(0) = 4 we obtain 



1 



1 



1 



1 



9 (1-w) 3 3 (1-w) 2 9 (1-w) 3 / 2 
4 1 35 



3 (1-w) 1 / 2 9 ' 



(B2) 



Expanding this generating function in a Taylor series 
then yields the result for (N 2 (N)) quoted in Eq. (f26|). 



APPENDIX C: GENERATING FUNCTION FOR 
FIRST MOMENT 

Here we solve the recursion formula Eq. ( |30|) for 
(Nk{N)}. We first introduce the generating function 
Xk{w) = Ew=i(^(^)) id 1-1 to eliminate the variable 
N and convert Eq. (^) into a differential equation that 
relates Xk and Xk-i- This equation is further simplified 
by making the transformation 



X k (w) = (1 - iu)* 1 Uk{u), u = 



1 



The resulting equation is 

^ = (fc-l)Z4-l, k>2. 
au 

Rewriting our previous solution ( p3| ) as 
2 



1. (CI) 



(C2) 



Ux{u) 



3 U 



2u l + 2u + 2, 



(C3) 



one can solve Eqs. (C2) subject to the initial condition 
U k (u = 0) = for k > 2. The final result is 



_ 4u fc+2 4u fe+1 2u k 

Uk{ - U > ~ k(k + i){k + 2) + kJk+T) + T 



2u 



fe-i 



Using the binomial formula, we transform Xk(z) into the 
series 

4 1 4 1 

Xk ^ W > ~ fc(fc + l)(fc + 2) (1-w) 2 + 3 (1-w) 1 / 2 

+ 2E(-ir^|( fc_1 ) (i--) (a - 1)/2 . 

4 a + 3 V. a / 

a— 1 x ' 

Expanding Xk(w) in a Taylor series in w we obtain 
(Nk{N)}. The analytic expressions for (Nf-(N)) with 
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k < 5 are obtained by expanding X k (w) in a Taylor se- 
ries. This gives 

2 4 T (N - -) 

1 4 r Civ — i) 
4 r (jv - |) 

" 50F T(N) ' 

12 rpv-j) 5 5 
r(A) -3^.1 + 3^ 

2 4 r Cat - -) 

(JV 5 (JV)) =-4iY +:r i^ -6(5 



105 J ' ' 30F r(iv) ""^ 
24 r (iV - I) 20 20 

50F r(iv) -y<W + T <W 
9 r(iv-|) 



7^ r(A r ) 

Generally, there are slightly different formulae for even 

k 

(N 2k (N)) = n 2k N + ^A k] 5 N3 

3 = 1 

1 4(i + i) /^fc - f\ r (AT - | - j) 



^ 2* + 3 V 2* y r(|-i) r(jv) 



and odd 



(N 2k+1 (N)) = n 2k+1 N + J2 B kj S Nj 
4(i + l) /2*\ T(7V 



+ & ^TT W r(|-i) r(jv) 

indices. Here the n k are given by Eqs. (||) and explicit ex- 
pressions for the coefficients A k j and B k j could be found 
by expanding the polynomials in the generating functions 
X 2k {w) and X 2k+1 (w). 



APPENDIX D: HIGHER MOMENTS 

Starting from Eq. (p9]), a straightforward computation 
yields 



(fc - ljjVfc-i + fciV fc 
2N 



(Dl) 



where the correlation function on the left-hand side is 
a function of N + 1 and those on the right-hand side 
are functions of N. Obviously, (N^) is coupled with 
(N k -iN k ). The recursion relation for this correlation 
function reads (for k > 3) 

(Nu-iN k ) =(l- 2 -^) <JV*_xJV*> + ^ (JVgLx) 



27V 



2iV 



(N k - 2 N k ) 



2N 



k-l 
2N 



k-l 



(D2) 



Fortunately no higher-order correlation functions ap- 
pear, and additionally the total index decreases, i.e., 
\N k ), whose total index is 2k, involves the correlation 
function (N k _iN k ), whose total index is 2k — 1. One 
therefore can determine all correlation functions by start- 
ing from the smallest total index and then working up to 
larger indices. For example, the first non-trivial corre- 
lation function (NiN 2 ) whose total index equals three 
satisfies an equa tion slightly different from the general 
form of Eq. (p2[), viz., 



(NiN 2 



+ 



2N 



2N 



(JV1JV2 



N 



(N 2 



(D3) 



Notice here that we already know (Nf). 

We can solve for (NiN 2 ) using the generating function 
technique. We define the generating function Zi(w) = 
Y,N>i( N i( N ) N 2(N)) w N ' x which satisfies the differen- 
tial equation 

2 il - w )^l = -Z 1 + y 1 +2w-^-, (D4) 
aw aw 



with solution 



1 



1 



1 



9 {l-wf 5 (l-w) 2 9 (1-iu) 3 / 2 
4 1 



47 /, M/2 35 

"8(13^-15 < 1 - W ^ + T- 
Expanding Z\ (w) in a power series in w we obtain 



(D5) 



» 1 , n 1 10 r jv + i 



4 r(JV-§) 47 T(N 



3v^F r(7V) 300? T(N) 
35 . 

Asymptotically, (NiN 2 ) — > (Ni)(N 2 ) as expected. 

There are two correlation functions, (A 7 !) and (NiN 3 ), 
whos e t otal index equals four. The former satisfies 
Eq. (pit) with fc = 2, i.e., 
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<^f> = (l-§) <^ 2 >+(^ 



+ 2N 2 + 2iViiV 2 
2N 



from which we determine the generating function 



18 (1 - wf 10 (1 - w) 2 9 (1 - w) 3 / 2 
4 1 94. , 1/2 49 55 

+ 9 (l-t P )Va -iS (1 " w) + T~iS V> ' 



Expanding 3^2 ( w ) we obtain 



(N^N)) = -N(N + l) + -N + 



4 r(iv + i) 



36 



4 r TV 



10 



+ 



9^ T(N) 



47 r (TV 



90F T(N) 15VtF r(/V) 
49 55 

In the large iV limit, we find that variance grows linearly 
with N according to a 2 ~ N. It appears that 



al -»• /i fc TV 



as TV — > oo, 



(D6) 



for all fc, although wc solved only the cases k = 1 and 2, 
l 

9 



where /ii = g and ^2 = j§q- 



APPENDIX E: SCALING FUNCTION IN THE 
DISCRETE APPROACH 

To extract the scaling function from the generating 
function J\f(w, z) we adapt the technique employed in 
Sec. IV for discrete variables. We first write 



1 + S y/l — W 



(El) 



and keep s finite while taking the w — > 1 limit. We again 
consider the modified generating function 



d_ 

Op: 



OO OO / 7 \ 

JV=lfc=l Wjv/ 



On the right-hand side of this equation we have already 
replaced (k + 2)(k + l)k(N k (N)) by 4NF(k/\/N) as im- 
plied by Eqs. (§)-(§). 

Substituting the exact expression (^) for the gener- 
ating function into the left-hand side of Eq. (E2) and 
keeping only the dominant contribution gives 



4(l-w)-^ 2 J{s), 



(E3) 



with J(s) given by Eq. (pL5[) . To simplify the right-hand 
side of Eq. (E2) we substitute Eq. (El) and replace the 
sums by integrals. The dominant contribution in the 
w — > 1 limit is 

/>oo p 00 

4(l-w)- 5/2 d£e-t s d^-qF^-^^e^ , (E4) 
Jo Jo 



where £ = k^/l — w and 77 = N(l — w). Therefore the 
double integral in Eq. (|E4| ) is equal to J(s). The dou- 
ble integral can be interpreted as the Laplace transform 
$(s) = J °° d£ exp(— s£)$(£) of the function 



*(£)=/ d^F^-^^e- 71 . (E5) 
Jo 

We already know how to solve 3>(s) = J(s), so 



<T> (i , = M -^i(l + y) - ; 



(E6) 



To determine -F(£), we must solve the integral equation 



(E6) with $(£) given by Eq. (|Eq). To solve this integral 
equation, notice that <&(£) is almost a Laplace transform 
of function F. Indeed, if instead of r\ and Fi^v^ 1 ! 2 ') we 
use C and G(£) defined according to 



C = |f, g(0 = cf(c 1/2 ), 



(E7) 



then we obtain $(£) = p 2 G(p), with p = £, 2 being the 
Laplace variable and G(p) = f£°d£G(£) exp(—p(). Re- 
writing the integral equation (Eq) in terms of p gives 



G{p) 



P 



P 



-3/2 



1 



1 



■ P 



;P 



-1/2 



2" 2 

Inverting this Laplace transform yields |2 



exp(-v/p) . 



G(C) = Cerfc 



1 



2C + 1 



-1/4C 



(E8) 



where erfc(a;) = ^7= dt exp(— t 2 ) is the complemen- 
tary error function. Since F(£) = £ 2 G(£~ 2 ), see Eq. (|E7|), 
we arrive at the scaled average degree distribution quoted 
in Eq. m. 
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